Stable isotopes unveil one millennium of domestic cat paleoecology in Europe

The domestic cat is the world's most popular pet and one of the most detrimental predators in terrestrial ecosystems. Effective protection of wildlife biodiversity demands detailed tracking of cat trophic ecology, and stable isotopes serve as a powerful proxy in dietary studies. However, a variable diet can make an isotopic pattern unreadable in opportunistic predators. To evaluate the usefulness of the isotopic method in cat ecology, we measured C and N isotope ratios in hundreds of archaeological cat bones. We determined trends in cat trophic paleoecology in northern Europe by exploiting population-scale patterns in animals from diverse locations. Our dataset shows a high variability of isotopic signals related to the socio-economic and/or geomorphological context. This points toward regularities in isotopic patterns across past cat populations. We provide a generalized guide to interpret the isotopic ecology of cats, emphasizing that regional isotopic baselines have a major impact on the isotopic signal.

www.nature.com/scientificreports/ from other sites, except natural sites (SI Appendix Table S9). This result was the effect of much lower δ 15 N values in specimens from rural (and natural) locations. Cats from seaports held a distinct position compared to other sites. These cats exhibited statistically significant different δ 15 N values to cats from other types of sites; except urban sites and castles. Their δ 15 N elevated values were attributed to higher contribution of marine food resources. For δ 13 C values, seaports statistically differed to strongholds, castles, and rural sites.
For the North Sea region, isotopic values significantly differed between most site types, in terms of both δ 13 C and δ 15 N (Figs. 4 and 5; SI Appendix Tables S10 and S11). In particular, there was a significant difference between urban sites and seaports (p < 0.001), which was not detected in the Baltic region (though, in the Baltic region, the null hypothesis of equal means was fulfilled for urban sites and seaports with weak probabilities, p < 0.15). The lowest δ 13 C and δ 15 N values were found at castle and rural sites, whereas the highest values related to seaports and a monastery located on the coast. The urban sites exhibited a wide range of values. The δ 15 N signal of single samples overlapped with the highest δ 15 N values detected at seaports, whereas the lowest signal fell in the range of castle and rural sites.
Cat isotope signal vs. geomorphological setting. To detect geomorphological factors that might shape the isotope signal of cats, we grouped the sites into five categories: inland, river valley, lakeshore, sea coast, and lagoon/fjord ( Fig. 4; Dataset S1; definitions presented in "Classification of the sites" section).
For cats from the Baltic region, we found no statistically significant difference between sites next to freshwater and brackish water bodies (lakes, lagoons, fjords, and larger rivers). This lack of differences was visible in terms of both δ 13 C and δ 15 N values (high p-values of the analysis of variance, usually p > 0.2, SI Appendix Tables S12 and S13). Inland and sea coast locations differed from each other and from all other geomorphological settings in terms of δ 13 C values. In contrast, δ 15 N values differed for certain pairings (e.g., inland vs. sea coast, inland vs. lagoon/fjord, sea coast vs. river valley, and sea coast vs. lagoon/fjord), with a weak similarity for inland vs. lakeshore pairs (p = 0.0535, SI Appendix Table S13). River valleys and inland settings had the broadest range of δ 15 N values (Fig. 5).
For the North Sea region, we detected statistically significant differences between most of site geomorphology categories in terms of both δ 13 C and δ 15 N values (analysis of variance p < 0.005, but was much lower in most pairs, SI Appendix Tables S14 and S15). The only statistically similar pair is sea coast vs. river valley, for δ 15 N values. However, p value for this pair is close to the 0.05 confidence level. Similar to the Baltic Sea region, cats from river valleys had the widest range of δ 15 N, as well as a wide range of δ 13 C values (Fig. 5d). . δ 13 C and δ 15 N values vs. specimen chronology in the cats analyzed in this study. The chronology reported for each specimen is the central point (± 1-year accuracy) of the site chronology range or the central point of the 95.4% probability range of the radiocarbon date. For full chronology ranges see Datasets S1 and S2 and for each socio-economic context shown separately see SI Appendix, Figs. S1 and S2.

Discussion
Potential sources of variability in the stable isotope values of cats. In general, the carbon and nitrogen isotopic signal (δ 13 C and δ 15 N) of a carnivore is controlled by many factors. These include biological factors (choice of prey, spectrum of prey, and isotopic variability of prey) and indirect environmental factors, which are responsible for the isotopic signal in prey (connected to vegetation type, aquatic/terrestrial ecosystem, soil development stage, climate, altitude) 26,27,[49][50][51][52][53] . For domestic carnivores, direct or indirect anthropogenic impacts might also play an important role in shaping isotopic signals. These include, but are not limited to, various human activities: (i) feeding on consumption/production waste or on selected food; (ii) creation of new landscapes, opening up habitats for new prey species; (iii) farming selected plants, which serve as a food base for herbivores (for pest rodents hunted by carnivores, or for domestic stock delivering meat and/or dairy that is further fed to carnivores); and (iv) supplying food that would otherwise be of limited availability to a terrestrial carnivore (such as food acquired from fishing or distant trade). The feeding behavior of cats is also driven by local conditions, including currently available food resources, their relationship with people (indoor, free-roaming, feral), and individual behavioral characteristics 1,4 . Recent studies of modern cat diet indicate high behavioral and feeding differentiation in domestic cats 1,17,18 . Terrestrial vs. aquatic resources. The most important factor shaping the isotope signal of the studied cats in the current study was the proportion of terrestrial and aquatic food resources in the diet. However, the regional isotope baseline (particularly the marine baseline), determined the degree to which differences could be detected through isotope values. In our dataset, the terrestrial and marine isotope signals were significantly different for cats from the North Sea region compared to those from the Baltic Sea region. Furthermore, the local background of human activity (especially cultivation of C 4 plants and location of the fishing grounds and the type of fish exploited) could enhance or obscure the isotope indicators of consumed resources. Cats from variable locations situated near water bodies (both marine and freshwater: sites on the sea coast, lagoons/fjords, river valleys, and lakeshores) had similar isotopic patterns with high δ 15 N values, whereas inland (i.e., far from the water bodies) cats had lower δ 15 N values (Figs. 5d and 6). Such elevated δ 15 N values are typical for aquatic ecosystems 53 . Thus, water-derived food resources contributed significantly to the diets of cats close to water bodies. Considering the socio-economic context of sites, higher δ 15 N values were expected in cats from seaports (Fig. 5b). Elevated δ 15 N values are also typical for cats from castles and strongholds in the Baltic region (Figs. 5b and 6). This is due to the cross-effect of cultural function and geomorphological setting of these sites, which were usually defensive locations situated on the edge of lakes or rivers (23 out of 25 of these sites in our database were attributed to river valley, lakeshore, or lagoon/fjord). www.nature.com/scientificreports/ Cats from non-coastal sites near to water (i.e., lakeshores, river valleys, and lagoons/fjords) were statistically dissimilar from those on the sea coast. This result was due to slightly lower mean δ 15 N values and clearly lower minimum δ 15 N values (Fig. 5d). This finding supported the expected difference between freshwater and marine fish 54,55 . In parallel, much lower δ 13 C values were obtained for cats from near-water non-coastal sites compared to cats from seashores (Fig. 5c). Despite being located quite far from the seashore, Bremen did not follow this relationship, as it had marine-like δ 13 C and δ 15 N values. However, availability of the marine resources in Bremen might have been due to the towns' trade hub function (see further discussion in "Cat isotope signal and the level of urbanization" and "Local conditions" sections).
The nitrogen and carbon isotopic signals in cats suggest an opportunistic shift toward local aquatic resources when available and a high reliance on such resources. This phenomenon might be related to scavenging fish remains or cats hunting terrestrial piscivores. Such prey might include sea birds, common in port towns and characterized by a marine-like δ 15 N signal 56 , and omnivorous rodents, which forage on fishery wastes and exhibit elevated δ 15 N values in port towns 57,58 . Cat isotope signal and the level of urbanization. Cats from urban sites and seaports had the most varied δ 13 C values (Figs. 5a and 6). Cats from seaports had the highest mean δ 15 N values and were characterized by the lowest variability in δ 15 N (Figs. 5b and 6). These differences in δ 13 C and δ 15 N values variability indicated different primary food sources. The nitrogen signal is mainly related to the input of terrestrial (low δ 15 N values) vs. aquatic (high δ 15 N values) resources. The carbon isotopic signal depends on marine (high δ 13 C values, but there are other differences between littoral and pelagic ecosystems) vs. non-marine (terrestrial and freshwater) resources. The presence of C 4 plants might also elevate the latter. The observed low δ 15 N and high δ 13 C variation in seaport cats indicated that aquatic resources were an important component of their diet, with mixed marine and non-marine resources, or mixed marine-littoral and marine-pelagic resources, or both.
The diet of cats from coastal locations was mainly driven by marine-derived food, whereas that of urban cats included a greater variety of food resources. Urban cats had high mean δ 15 N, but also, highly varied δ 15 N values (Fig. 5b). However, high-δ 15 N values of freshwater and/or marine food was still available to urban cats,  www.nature.com/scientificreports/ regardless of whether it was naturally present in the local ecosystem. This finding demonstrates the importance of anthropogenic resources for these cats, as the availability of water-derived food was certainly due to trade. Medieval logistics allowed transport of fresh sea-fish up to 150 km inland 59 . This heterogeneity in the diet of cats in urban areas has also been identified in studies of modern cats 11,32,58,60,61 : cats consume more anthropogenic food in urban habitats, while predation rate and prey diversity are highest in rural areas. Cats connected to natural locations exhibited relatively low variability for both carbon and nitrogen (Fig. 5a). They were characterized by low δ 15 N and δ 13 C values, which is typical for forest animals 10,34 . This observation indicates a narrow isotopic niche and trophic specialization, and shows that cats probably mainly inhabited forests and exploited wild prey. Studies of recent cats also support this explanation, showing that feral domestic cats tend to have a narrower trophic niche than house cats, with their diet similar to that of wildcats 18 .
The carbon and nitrogen isotopic ranges of rural cats overlapped with those from natural sites (Fig. 5a,b). However, rural cats exhibited a broader range of values, especially for δ 13 C. Like cats from natural sites, rural cats had lower δ 13 C values, and the lowest mean δ 15 N values, compared to other locations. These results indicate the use of just terrestrial prey resources, and predation over wild prey species. Increased predation rates in rural areas commonly occurs in modern cats 11,34,60,61 . While urban cats are more likely to be locked inside or living close to households, rural cats have lower population densities but have more extensive home ranges and live closer to natural habitats 11,17,18,62,63 . Local conditions. Detailed analysis of our dataset allowed us to detect some abnormalities in the isotopic pattern. For instance, several sites from the Baltic region (Kołobrzeg, Puck, Napole, and Visby) stood out from the whole dataset due to their high δ 13 C signal ( Fig. 4; Datasets S1 and S2). Several potential factors could be responsible for such a high δ 13 C signal: low forestation 26 , high elevation 64,65 , or significant contribution of marine food resources 54,55 . Another factor that elevates δ 13 C values might be attributed to the important contribution of C 4 plants (e.g., millet) at the basis of the trophic web. These C 4 plants were possibly stored cereals, which served as a source of food for pest rodents, which, in turn, were preyed upon by cats. The presence of millet was confirmed in archaeobotanical material from Puck 66 . However, in the case of coastal sites, such as Puck, Kołobrzeg and Visby, the marine resources could also play an important role. Cats had significantly different δ 15 N values in these seaports versus Napole. These seaports had the highest, whereas Napole had the lowest values out of the entire dataset. This difference followed the seaport vs. rural site function, and likely reflected differences in the proportion of marine food available at each site.
Another interesting example concerns fortified sites, namely strongholds and castles from the Baltic region. Cats from these sites exhibited the lowest isotopic variability for δ 13 C values and moderate variability for δ 15 N values. The absolute δ 13 C values were similar to natural and rural sites, whereas δ 15 N values were usually higher (Fig. 4). This result indicates that freshwater fish formed an important component of the diet, which is logical considering the location of these sites (23 out of 25 sites were attributed to the river valley, lakeshore, or lagoon/ fjord geomorphological setting). Cats inside these closed settlements might have preferred living close to humans, and therefore had an easy access to table refuse.
The castle of Logne was the only castle situated in the North Sea region in our database. Logne cats had a clearly lower δ 15 N signal compared to cats from castles in the Baltic Sea region (Fig. 4). This result was unexpected when considering the generally higher δ 15 N values baseline in the North Sea region. Cats from this location (which is rather far from the coast), might have had limited access to fish and other water resources. The region has been densely forested until recent times, which suggests that the Logne cats could have had easy access to forest resources (via hunting or access to the table refuse of game animal meat), reflected in their low δ 15 N and δ 13 C values.
The opposite situation is visible in cats from Bremen, recording marine-like isotopic signals despite its location in the freshwater zone of the Weser 67 . Marine fish, however, were commonly available in Bremen as evidenced by zooarchaeological studies and historical records. Some anadromous species (e.g., sturgeon, salmon, sea trout, smelt, twaite shad, and flounder) used to enter the Weser during their spawning season and might have been fished locally 68 . The Bremen fishermen were organized in a guild from at least 1489, which held the fishing rights down to the river mouth 69 . Additionally, Bremen was an important trade hub since at least the eleventh century, particularly involved in the stockfish and herring trade 70 .
A model of cat isotopic ecology. In general, domestic cats may follow three main foraging strategies: (i) hunting, (ii) scavenging, mainly on human waste, and iii) direct feeding by humans. Considering these strategies, we may define the following factors responsible for an isotopic signal of cat diet, synthesized in Fig. 7: • Presence of cereal crops, cultivated in open fields, characterized by elevated δ 13 C values 26,27 , and available to synanthropic rodents (such as mice and rats, hunted by cats) either directly from fields, storages in barns or granaries, or from processed products; • Cultivation or import of C 4 cereal crops with high δ 13 C signals (such as millet or maize) 50 www.nature.com/scientificreports/ • Access to fish (and other marine resources), usually characterized by an elevated δ 15 N signal and variable δ 13 C values, depending on fish ecology and physiology (i.e., differences between freshwater and marine fish 54,55 , geography (e.g., the difference between the Baltic Sea and North Sea 47,48 , or between pelagic and littoral waters 53 ).
In a human-dominated landscape, the isotopic signal of cat diet can be impacted by all of the listed factors. However, certain basic patterns related to variable anthropogenic and geographical settings can be predicted. In rural areas, crop pests are easy prey; however, the proximity of natural environments encourages the opportunity for cats to wander and feed on wild prey. Modern rural cats have a more diverse diet and a higher predation rate compared to those inhabiting urban sites 11,60,61 . Towns favor larger populations of synanthropic rodents, due to food stores and refuse. In addition, urban areas as trading spots provide various imported food resources, including marine fish, even in inland locations. However, even the largest inland towns cannot compete with seaports in terms of the availability of marine resources.

Conclusions
The ecology of the domestic cat is complex because it is an opportunistic exploiter that can maneuver among multiple available ecological niches 1 . Modern cats sharing the same environment, especially town-dwellers, exploit highly variable habitats and, hence, have highly varied trophic ecology 1,12,16,17,31,72 . Some cats may live as "backyard cats," (i.e., welcome but free-roaming mouse hunters), which are occasionally fattened by their human neighbors. Others might be "feral outcasts," which are the barely tolerated semi-wild dwellers of cellars, roofs, and hovels, living on their own and avoiding human contact. Still others might live as "indoor pets." Some cats even, exploit all these lifestyles.
The individualistic behavior and the taphonomic limitations mean that the isotopic signature of a single fossil cat specimen, or even of several specimens, might not represent the entire population. These limitations were revealed in our study through the high variability of the isotopic signal within the cats from some socioeconomic types of sites (particularly those from seaports, urban and rural settlements) and geomorphological settings (especially river valleys and sea coasts). However, our data also revealed relatively low isotopic variability in cats from specific locations (mainly natural sites, outside of cultural context). This phenomenon most likely relates to feral cats, which can permanently or temporarily adopt a wildcat-like lifestyle. Similarly to modern cat populations, the medieval cats tended to have more heterogenous diet in urban areas, as they consumed more anthropogenic food. At the same time, the predation rate was higher in rural and natural areas.
The observed isotopic diversity means that stable isotope analysis could be a useful tool in elucidating the diet and lifestyle of individual cats. However, a local/regional isotopic baseline must always be carefully considered, as no global range may be arbitrarily attributed to the particular ecology of a cat. This study shows a promising perspective concerning the trophic ecology of fossil and recent cats, offered by an analysis of diet-related stable www.nature.com/scientificreports/ isotopes preserved in bones. Further investigations using sulfur stable isotope analysis and compound-specific carbon and nitrogen isotope analysis may deliver more in-depth insight into the domestic cat ecology.

Materials and methods
Zooarchaeological material. We collected the bones of small felids from archaeological and paleontological sites located in Poland, Germany, and Belgium ( Fig. 1; Datasets S1 and S2). Our dataset included various geographical locations, covering human settlements of various socio-economic contexts, and natural sites. The available remains were either single bone specimens or almost complete skeletons. In total, 190 bones originating from 60 archaeological sites at 40 locations were sampled. The characteristics of the sampled specimens and sites included in this study are listed in (Datasets S1 and S2). We also included literature data of 18 felid bones from Sweden, Norway, UK, Belgium and Germany, from previously published studies [38][39][40][41][42][43][44][45] (Dataset S3). Bones most probably belonging to the same individual (based on: the same stratigraphic unit; different skeletal elements; similar size, individual age and preservation state), were sampled and analyzed separately, but their isotopic values were averaged (Dataset S4).
Taxonomic identification-domestic cat vs. European wildcat. Preliminary identification to the genus level (Felis sp. Linnaeus, 1758) was carried out in the Royal Belgian Institute of Natural Sciences and the Institute of Archaeology, Nicolaus Copernicus University in Toruń, using reference collections of domestic cat and European wildcat skeletons. The bones were measured according to the methodology by A. von den Driesch 75 , using a digital caliper (0.1 mm reading error) (Dataset S2). To determine the species (domestic cat, F. catus Linnaeus, 1758 vs. wildcat, F. silvestris Schreber, 1777, considered here to be European wildcat, F. silvestris silvestris) we adopted the discrimination method based on relative size using the log-ratio technique by T. P. O'Connor 76 (Datasets S5, S6 and S7). The reference standards for domestic cat and European wildcat size ranges were based on Z. Kratochvíl's database 77,78 .
Due to limitations in distinguishing domestic cat bones from the European wildcat based on osteometric and morphological features, we also used mtDNA analysis, especially for specimens where the size ranges of domestic cats and European wildcats overlapped (SI Appendix Fig. S3; Dataset S2).
On the basis of morphometry and/or mtDNA analysis, we distinguished 139 domestic cats and 12 European wildcats. The specimens for which we could not obtain genetic or morphometric taxonomical information (n = 31) likely represent domestic cats due to their relatively small size, and have been regarded as such in further analyses. Specimens extracted from previous publications (n = 18) were also considered domestic cats.
Cat age criteria. We selected bones of individuals that were older than the nursing period, to avoid the milk diet obscuring trophic ecology of individuals. Cat dependence on the mother's milk is short and weaning is considered finished by about 7 weeks after birth 1 . Thus, for our study we selected adult individuals (with bone fusion completed) and subadult individuals, i.e. older than weaning age, based on combined bone fusion or tooth eruption 79 . For fragmented bones, when it was not possible to apply the criteria above, we used overall bone size and stage of bone ossification as indicators of post-weaning age (Dataset S2).
Classification of the sites. Chronology. The chronology of the specimens was presented as ranges in years AD. These ranges were based on the radiocarbon determination of the specimen if available, or/and on stratigraphic and archaeological dating (Datasets S1 and S2), using the published data (and in some cases, unpublished archaeological reports) (Dataset S2).
We used all radiocarbon dates available in the published literature for the previously studied specimens 8,10,80 , and new data obtained during the course of current projects by the authors (Dataset S2). In all cases, the dated fraction used was bone collagen. Most dates were obtained in the Poznań Radiocarbon Laboratory (Poznań, Poland), while some dates were obtained in the Royal Institute for Cultural Heritage (RICH, Brussels, Belgium). The presented chronology range was the end points of the 95.4% probability range, achieved through calibration in OxCal, v. 4.4 81,82 . The ranges were provided with a ± 1-year accuracy.
The resolution varied among sites, from intervals of a part of a century to several centuries-long. We attributed chronological ranges to the stratigraphic data provided as centuries or parts of centuries in the following way: • a century: AD 1X00-1Y00; • early century: AD 1X00-1X25; • late century: AD 1X75-1Y00; • mid-century: AD 1X50; • 1st half of a century: AD 1X00-1X50; • 2nd half of a century: AD 1X50-1Y00.
On rare occasions the chronological data were more precise (i.e., narrowed to a period before or after a particular year). In such cases, we used these years as the range end. In cases where chronology was provided in a descriptive way, we adopted the following scheme: late medieval = AD 1250-1500; post-medieval = AD 1500-1800.
Geomorphology. We attributed a geomorphological context to the studied sites based on their geographic location (Dataset S1). We used the following types of geomorphological setting: www.nature.com/scientificreports/ • sea coast (less than 5 km from the coast of a sea or a bay; taking into account home size ranges and mobility to acquire food 1,83,84 ); • lagoon/fjord (over 5 km from an open sea coast, but less than 5 km from a brackish basin, such as a lagoon or fjord, which was weakly connected with the open sea); • lakeshore (less than 5 km from a larger lake, or at an island); • river valley (within a valley or at the edge of a valley of a large river); • inland (none of the above, i.e., over 5 km from the nearest sea coast, brackish basin, lake or large river valley).
Site socio-economic context. We attributed a socio-economic context to each site based on data from the published literature (Dataset S1), also taking into account the chronology of the studied specimens for sites where the function changed during medieval/post-medieval times. We attributed the following functions: • natural (no historical nor archaeological indications for any anthropogenic use of the site during the considered period; mostly cave sites); • rural (small settlement, possibly farmland, often loosely constructed with wooden buildings); • urban (towns, administrative centers, densely built-up, densely populated); • seaport (important marine ports; all of them were also urban sites); in this category, we also include Bremen, which was involved in marine fishing and trade in the past 68-70 ; • castle (late medieval fortified masonry seat of a noble family, a bishop, or knights of the Teutonic Order); • stronghold (mostly early medieval fortified settlement or refuge with wooden and earthen embankments; only in the eastern part of the studied area); • monastery (late medieval ecclesiastic settlement; only one example, Koksijde Monastery in Belgium).
Latitude, longitude, and distance from the sea. For cave sites, we provided the geographic coordinates of cave entrances. For open-air sites, we provided the coordinates of the center of a village/small town, center of a castle, or center of a street. We used Google Maps to obtain coordinates (decimal N and E values, ± 0.00001°). To provide the distance from the sea, we used the "measure distance" option in Google Maps. We measured this as the shortest distance from the coordinate-collecting point to the nearest sea coast line (± 1 km).
Stable isotope analysis. The bulk stable isotope analysis of carbon and nitrogen was performed on collagen extracted from bones. Dentine was not sampled to avoid bias from inter-tissue fractionation.
For samples with lab. no. CAT, collagen purification method followed the protocol after Krajcarz et al. and Bocherens et al. 10,85 . Isotopic measurements (δ 13 C and δ 15 N) and elemental measurements (elemental C:N ratio, calculated into atomic C:N ratio) were performed in duplicate at the Stable Isotopes Laboratory at the Institute of Geological Sciences, Polish Academy of Sciences (Warszawa, Poland) using a Flash EA 1112HT elemental analyzer (Thermo Scientific) connected to a Delta V Advantage mass spectrometer (Thermo Scientific). Mean SEs were < 0.33‰ for δ 13 C and < 0.43‰ for δ 15 N values. We normalized all measurements to δ 13 C values of USGS40 and USGS41 standards and all δ 15 N values to the IAEA 600 standard. Details are provided in SI Appendix Supplementary Text and Table S16.
For all of the remaining samples, collagen extraction was carried out at the Royal Institute for Cultural Heritage (RICH) in Brussels. Samples were prepared according to the method of Longin 86 , which was modified with an additional NaOH-wash between demineralization and hydrolyzation 87,88 . Collagen subsamples were sent to the Division of Soil and Water Management of the KU Leuven for stable isotope (δ 13 C and δ 15 N) and elemental (%C, %N) analyses. The atomic C:N ratio was then calculated. Analyses were carried out on a Thermo Flash HT/ EA elemental analyzer, linked to a Thermo DeltaV Advantage Isotope Ratio Mass Spectrometer via a ConfloIV interface (Thermo Scientific). Data calibrations were done using IAEA-600 (caffeine) and two in-house standards (Leucine and muscle tissue of Pacific Tuna), which were previously calibrated versus certified standards. Analytical precision was typically better than 0.15‰ for both δ 13 C and δ 15 N values, this was based on the average of four replicates of the different standards in each batch of samples.. Details are provided in SI Appendix Supplementary Text and Table S16.
Data were reported in delta (δ) notation: δ 13 C or δ 15 N = ([R sample /R international reference ] -1) × 1000 in ‰ (parts per thousand), where R is 13 C/ 12 C or 15 N/ 14 N. The international references were Vienna Pee Dee Belemnite (VPDB) for carbon and atmospheric nitrogen (AIR) for nitrogen.
Statistics. We used one-way Analysis of Variance (ANOVA; or alternatively Kruskal-Wallis or Welch F tests in the case of non-parametric distributions or unequal variances) to check for statistical differences between isotopic variances of selected groups of samples (such as samples representing: chronological intervals, site cultural functions or site geomorphologies). The homogeneity of variance and normality were tested with Levene's tests (from means and medians) and a set of normality tests (Shapiro-Wilk, Anderson-Darling, Lilliefors and Jarque-Bera tests) 89 . Because in many cases the groups exhibited non-normal distribution and unequal variances, post-hoc pairwise comparisons were difficult to perform. Therefore, analysis of variance was run separately for each pair. In cases where post-hoc tests were possible for the entire set of pairs, we applied Tukey's test (after a significant ANOVA) with the Bonferroni correction of p values for multiple comparisons. All tests were performed with PAST software, v. 4.05 90 . www.nature.com/scientificreports/ working with ancient DNA following strict rules to avoid contamination. Samples were prepared and processed as described by Baca et al. 8 . DNA was extracted using silica-coated magnetic beads according to Rohland et al. 91 . Double-strand, double-indexed library preparation, target enrichment of mitochondrial DNA (mtDNA), sequencing on the Illumina platform, and processing of the sequencing reads were performed as in Baca et al. 8 , with the modification that only positions with a minimum of 3× coverage were called. We obtained mtDNA sequences for 57 specimens, which were used in phylogenetic analysis, along with 123 already published sequences 5,8,10 . The Maximum Likelihood Tree was generated using IQtree 92 . The implemented ModelFinder determined that GTR is the best-fit model and the branch support was assessed using ultrafast bootstrap method (SI Appendix Fig. S3). All mtDNA sequences obtained in this study were deposed in GenBank under accession numbers: OL654306-OL654362.

Data availability
All data are available in the main text or the supplementary materials. DNA sequences obtained in this study were deposited in GenBank under accession numbers: OL654306-OL654362; direct link: https:// www. ncbi. nlm. nih. gov/ popse t/? term= 22599 73155. www.nature.com/scientificreports/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.